  SUBROUTINE INIT_VAR &
  ( M, Mice, ns, ms, ml, Tinitlength, init_T, skin, zero_model, &
  & h10, l10, hs10, ls10, tempair, Ts0, Tb0, Tm, Tbb0, &
  & Sals0, Salb0, us0, vs0, h_ML0, &
  & rosoil, por, &
  & zTinitprof, Tinitprof, &
  
  & flag_snow, flag_snow_init, itop, nstep, &
  & h1, l1, hs1, ls1, &
  & veg, snmelt, cdmw2, velfrict_prev, &
  & roughness, eflux0_kinem, Elatent, &
  & totalevap, totalmelt, totalprecip, totalwat, totalpen, &
  & time, dhwfsoil, dhw, dhw0, dhi, dhi0, dls0, &
  
  & E1, eps1, Tsoil1, wi1, wl1, Sals1, &
  & rootss, qsoil, TgrAnn, qwater, &
  & oxyg, carbdi, Sal1, u1, v1, Tw1, Tskin, T_0dim, &
  & Ti1, Tis1, z_full, ddz, &
  & dz, T, wl, dens, lamw, &
  & dzeta_int, zsoil, pressure)
  
  ! Subroutine INIT_VAR initializes the prognostic variables of the model
  
  use PHYS_FUNC, only : &
  & MELTPNT, &
  & UNFRWAT, &
  & WL_MAX, &
  & WI_MAX, &
  & HENRY_CONST
  
  use PHYS_CONSTANTS, only : &
  & row0, g, Kelvin0
  
  use METH_OXYG_CONSTANTS, only : &
  & ch4_atm0, o2_atm0, co2_atm0, &
  & rel_conc_ebul_crit, &
  & Henry_const0_ch4, &
  & Henry_temp_dep_ch4, &
  & Henry_temp_ref
  
!  use ARRAYS, only : &
!  & qwater2, qsoil2 ! two-meth
  
  implicit none
  
  real(8), parameter :: E_init = 10.d0**(-5.5)
  real(8), parameter :: eps_init = 1.d-9
  
  ! Input variables
  integer, intent(in) :: M, Mice, ns, ms, ml, Tinitlength ! Grid dimensions
  integer, intent(in) :: init_T, skin, zero_model ! Driving parameters
  
  real(8), intent(in) :: h10, l10, hs10, ls10
  real(8), intent(in) :: tempair, Ts0, Tb0, Tm, Tbb0
  real(8), intent(in) :: Sals0, Salb0
  real(8), intent(in) :: us0, vs0
  real(8), intent(in) :: h_ML0
  real(8), intent(in) :: rosoil(1:ns), por(1:ns)
  real(8), intent(in) :: zTinitprof(1:Tinitlength), Tinitprof(1:Tinitlength)
  real(8), intent(in) :: zsoil(1:ns+2)
  real(8), intent(in) :: pressure
  real(8), intent(in) :: dzeta_int(1:M+1)
  real(8), intent(in) :: ddz(1:M)
  
  ! Output variables
  integer(4), intent(out) :: flag_snow, flag_snow_init, itop
  integer(4), intent(out) :: nstep
  
  real(8), intent(out) :: h1, l1, hs1, ls1
  real(8), intent(out) :: veg
  real(8), intent(out) :: snmelt
  real(8), intent(out) :: cdmw2, velfrict_prev
  real(8), intent(out) :: roughness
  real(8), intent(out) :: eflux0_kinem, Elatent
  real(8), intent(out) :: totalevap, totalmelt, totalprecip, totalwat, totalpen
  real(8), intent(out) :: time
  real(8), intent(out) :: dhwfsoil
  real(8), intent(out) :: dhw, dhw0, dhi, dhi0, dls0  
  
  real(8), intent(out) :: E1(1:M+1), eps1(1:M+1)
  real(8), intent(out) :: Tsoil1(1:ns), wi1(1:ns), wl1(1:ns), Sals1(1:ns)
  real(8), intent(out) :: rootss(1:ns), qsoil(1:ns), TgrAnn(1:ns), qwater(1:M+1)
  real(8), intent(out) :: oxyg(1:M+1), carbdi(1:M+1)
  real(8), intent(out) :: Sal1(1:M+1)
  real(8), intent(out) :: u1(1:M+1), v1(1:M+1)
  real(8), intent(out) :: Tw1(1:M+1), Tskin(1:2)
  real(8), intent(out) :: Ti1(1:Mice+1), Tis1(1:Mice+1)
  real(8), intent(out) :: z_full(1:M+1)
  real(8), intent(out) :: dz(1:ms), T(1:ml), wl(1:ml), dens(1:ms)
  real(8), intent(out) :: lamw(1:M)
  real(8), intent(out) :: T_0dim  
  
  ! Local variables
  real(8), parameter :: Ct_d_Ch = 2. ! Using results by West & Plug, 2007, may be estimated
                                     ! from 1.7 to 2.8
  
  real(8) :: z, h_ML0zv, h_talik
  real(8) :: Ti10, Ti11
  integer(4) :: i, i_ML
  integer(4) :: n_1cm, n_5cm
  logical :: flag
      
  h1 = h10
  l1 = l10
  hs1 = hs10
  ls1 = ls10
  do i = 1, M 
    E1(i) = E_init !+ float(i)/float(M)*(1.d-10  -  1.d-8)
    eps1(i)= eps_init !+ float(i)/float(M)*(1.d-18  -  1.d-14)
  enddo

  
  if (h10 > 0. .and. Tb0 > 0. .and. Tbb0 < 0.) then
  
!   Initializing the talik depth, assuming thermokarst lake, following West and Plug, 2007
    h_talik = Ct_d_Ch*h10
    
    do i = 1, ns
      if (zsoil(i) >= h_talik) then
        i_ML = i ! i_ML is the number of the first level below the estimated talik depth
        exit
      endif
    enddo
    
!   Temperature distribution in the talik
    do i = 1, i_ML
      Tsoil1(i) = Tb0 + float(i-1)/float(i_ML-1)*(0.d0 - Tb0)
    enddo

!   Temperature distribution below the talik
    do i = i_ML+1, ns
      Tsoil1(i) = 0.d0 + float(i-i_ML)/float(ns-i_ML)*(Tbb0 - 0.d0)    
    enddo

  else
  
!   Linear temperature profile in soil
    do i = 1, ns
      Tsoil1(i) = Tb0 + float(i-1)/float(ns-1)*(Tbb0-Tb0)
    enddo 
    
  endif


  do i = 1, ns
    Sals1(i) = Salb0 ! assuming, that salinity in ground is the same 
                     ! as one in near bottom layer of water
    if (Tsoil1(i) > MELTPNT(Sals1(i))) then
      wi1(i) = 0.
      wl1(i) = WL_MAX(por(i),rosoil(i),0.d0)-0.01
    else
      wl1(i) = UNFRWAT(Tsoil1(i),i)
      wi1(i) = WI_MAX(por(i),rosoil(i))-wl1(i)-0.01
    endif
  enddo
  
! Methane variables initialization  
  rootss(1:ns) = 0.d0 ! No roots
  veg = 0.d0 ! No vegetation at the bottom
  qwater(1:M+1) = ch4_atm0 ! Atmospheric concentration
!  qwater2(1:M+1,1,1:2) = 0.5*ch4_atm0 ! two-meth
  TgrAnn(1:ns) = Tsoil1(1:ns)
  
  qsoil(1) = ch4_atm0
!  qsoil2(1,1:2) = 0.5*ch4_atm0 ! two-meth
  do i = 2, ns
    qsoil(i) = rel_conc_ebul_crit*por(i)*(pressure + row0*g*h10) * &
    & HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, &
    & Henry_temp_ref, Tsoil1(i) + Kelvin0)
!    qsoil2(i,1:2) = 0.5*qsoil(i) ! two-meth
  enddo

! Carbon dioxide initialization
  carbdi(1:M+1) = co2_atm0 ! Atmospheric concentration

! Oxygen initialization  
  i_ML = int(M*h_ML0/h1)
  do i = 1, M+1
    if (i <= i_ML) then
      oxyg(i) = o2_atm0 ! Atmospheric concentration
    else
      oxyg(i) = o2_atm0 - o2_atm0*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
    endif
  enddo
  
  Sal1(1:max(i_ML,1)) = Sals0
  u1(1:max(i_ML,1)) = us0
  v1(1:max(i_ML,1)) = vs0
  do i = max(i_ML+1,2), M+1
    Sal1(i) = Sals0 + (Salb0-Sals0)*float(i-i_ML)/float(M+1-i_ML)
    u1(i) = us0 + (0.d0 - us0)*float(i-i_ML)/float(M+1-i_ML)
    v1(i) = vs0 + (0.d0 - vs0)*float(i-i_ML)/float(M+1-i_ML)
  enddo
 
  if (init_T == 2) then
    h_ML0zv = (Tm-0.5*(Tb0+Ts0))/(Ts0/h1-0.5*(Ts0+Tb0)/h1)
    i_ML = int(M*h_ML0zv/h1)  
  endif 
  if (init_T == 1 .or. init_T == 2) then
    Tw1(1:max(i_ML,1)) = Ts0
    do i=max(i_ML+1,2), M+1 
      Tw1(i) = Ts0 + (Tb0-Ts0)*float(i-i_ML)/float(M+1-i_ML)
    enddo
! The initial temperature profile is given from the input file
  elseif (init_T == 3) then
    do i = 1, M+1
      z_full  (i) = dzeta_int(i)*h1
    enddo
    call LININTERPOL (zTinitprof,Tinitprof,Tinitlength, &
    & z_full,Tw1,M+1,flag)
    if (.not.flag) then
      print*, 'The error while interpolating the initial &
      &temperature profile: terminating program'
      STOP
    endif
  endif
  
  if (skin == 0) then
    Tskin(1:2) = 0.d0
  else
    Tskin(1) = Tw1(1)
  endif  
  
! Initial temperature for zero-dimensional model 
  if (zero_model == 1) then
    T_0dim = Tw1(1)*0.5*ddz(1) + Tw1(M+1)*0.5*ddz(M)
    do i = 2, M
      T_0dim = T_0dim + Tw1(i)*0.5*(ddz(i-1) + ddz(i) )
    enddo
  else
    T_0dim = 0.
  endif
 
  Ti1 = 0.
  Tis1 = 0. 
  Ti10 = dmin1(tempair,-1.d-1) ! Initial ice surface temperature, Celsius
  if (l1 /= 0) then
    if (h1 > 0.) then
!     Ice over water 
      Ti11 = MELTPNT(Sals0)
    else
!     Ice over soil
      Ti11 = Tb0
    endif
    do i = 1, Mice+1
      Ti1(i)= Ti10 + (Ti11 - Ti10)*float(i-1)/float(Mice) ! Linear profile, if dzetai-grid is regular,
    enddo                                                 ! water layer underneath is assumed to exist
  endif
  
  if (hs1 == 0.) then   
    flag_snow = 0
    flag_snow_init = 1
    itop = 1
  else
    flag_snow = 1
    flag_snow_init = 0
    hs1 = dmax1(2.d-2, hs1)
    hs1 = int(hs1/0.01)*0.01
    n_5cm = int(hs1/0.05)
    n_1cm = int((hs1-n_5cm*0.05)/0.01)
    if (n_1cm == 0) then
      n_1cm = 1 
      hs1 = hs1 + 0.01
    endif
    itop = ms - (n_1cm + n_5cm)
    do i = itop, itop + n_1cm - 1
      dz(i) = 0.01
      T(i) = dmin1(tempair,-5.d0)
      wl(i) = 0
      dens(i) = 150.
    enddo
    do i = itop + n_1cm, ms-1
      dz(i) = 0.05
      T(i) = dmin1(tempair,-5.d0)
      wl(i) = 0
      dens(i) = 150.
    enddo
    dz(ms) = 0.5d0*dz(ms-1)
    T(ms) = dmin1(tempair,-5.d0)
    wl(ms) = 0.d0
    dens(ms) = 150.d0
  endif
 
  snmelt = 0.d0 
  cdmw2 = 1.d-15
  velfrict_prev = 1.d-2 
  roughness = 1.d-3
  eflux0_kinem = 0.d0
  Elatent = 0.

  totalevap = 0 
  totalmelt = 0 
  totalprecip = 0. 
  totalwat = 0
  totalpen = 0
  time = 0
  nstep = 0
  dhwfsoil = 0.
  dhw = 0.
  dhw0 = 0.
  dhi = 0.
  dhi0 = 0.
  dls0 = 0.

  lamw(:) = 1.d+3
  
  END SUBROUTINE INIT_VAR